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ABSTRACT 

Aims. We describe a new methodology to analyze the reionization process in numerical simulations: the chronology and 
the geometry of reionization is investigated by means of merger histories of individual HII regions. 
Methods. From the merger tree of ionized patches, one can track the individual evolution of the regions properties such 
as e.g. their size, or the intensity of the percolation process by looking at the formation rate, the frequency of mergers 
and the number of individual HII regions involved in the mergers. We apply the merger tree technique to simulations 
of reionization with three different kinds of ionizing source models and two resolutions. Two of them use star particles 
as ionizing sources. In this case we confront two emissivity evolutions for the sources in order to reach the reionization 
at z ~ 6. As an alternative we built a semi-analytical model where the dark matter halos extracted from the density 
fields are assumed as ionizing sources. 

Results. We then show how this methodology is a good candidate to quantify the impact of the adopted star formation 
on the history of the observed reionization. The semi-analytical model shows a homogeneous reionization history with 
'local' hierarchical growth steps for individual HII regions. On the other hand auto-consistent models for star formation 
tend to present fewer regions with a dominant region in size which governs the fusion process early in the reionization at 
the expense of the 'local' reionizations. The differences are attenuated when the resolution of the simulation is increased. 
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1. Introduction 

The first generation of stars appears after the period of dark 
ages (between z ~ 1090 and z ~ 20), creating a multitude 
of distinct HII regions in the Universe. These regions even- 
tually merge with other ionized regions, leading to a total 
reionization of the Universe. The diffusion of CMB pho- 
tons on the electrons released during reionization (see e.g. 



Komatsu et al. 2009) and the absorptio n features in the 
spe ctra of high red shift quasars (see e.g. Fan et al. 2006 
[Willott et al.|2007 ) suggest that the end of reionization oc- 
curred at z ^ 11 — 6 , with an average neut ral fractio n of 
hydrogen between 10 ~^ and 10~^ (see^Fan et al.||2006 ). 
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distribution and ionization fields are derived from the den- 
sity fields. 

In this context, one challenge is to describe the geometry 
and the time sequence of the reionization as they will be 
influenced by e.g. the formation rate of ionizing sources, 
the distribution of their formation sites or the size and 
growth of HII regions. Aiming at this description, many 
theoretical studies have thus been conducted focusing on 
the global evolution of physical fields such as the ionized 



fraction or temperatu re (see e.g. Aubert fc Teyssier 2010 
Finlator et al. 2009b). Alternately, many groups have un- 



ing on the properties of individual HII re^ 


^ions (Iliev et al. 


2006b| jZahn et al. 2007| [McQuinn et al. 


2007, 


Lee et al. 


2008 


Croft fc Altay||2008| 


Shin et al.[|200g 


, Friec 


rich et al. 


201ip. 



With the current article, we also aim at investigating 
the chronology and geometry of the reionization through 
numerical simulations. In particular, we describe the pro- 
cess through the study of the merger process of HII regions 
and for this purpose we track the individual histories of the 
ionized regions formed in numerical simulations. Using this 
point-of-view, we combine the analyzis in terms of individ- 
ual HII regions with the temporal evolution of global fields. 
As a consequence, it results in a 'local' perspective with 
a set of histories of reionization^ allowing to study their 
scatter throughout a cosmological volume. This will be an 
alternative approach to those already undertaken. 
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In this paper, we first focus on the procedure to sep- 
arate the different HII regions and to track each of them 
along tim e. For this purpose, we follow man y ot her inves- 
tigations dlliev et al.| |2QQ6bl |Shin et al.pQQSl and [Friedrich 
etaI]|2Qlip by using a Friends- of -friends algorithm (FOF) 
to characterize the different ionized regions. Then, to follow 
the HII regions' properties with time or redshift, we build 
their merger tree. Using the merger tree and investigat- 
ing its properties to infer the redshift evolution of the HII 
regions' formation number, their size evolution and their 
merger history, we suggest a new way to constrain the evo- 
lution of the reionization process. 

We test the method on three models for the ionizing 
sources in order to quantify their impact on the simulated 
reionization history. We also quantify the impact of the spa- 
tial resolution of the simulations on the history by perform- 
ing simulations with sizes of 200 and 50 Mpc/h boxes. By 
comparing the resulting reionization histories, we show how 
different source prescriptions lead to different properties on 
the evolution of HII regions as reionization progresses. 

This paper is organized as follows. In Section [2j we 
present the tools developed to investigate the history of 
the reionization in simulations such as the FOF algorithm 
and the merger tree of HII regions. In Section |3) we detail 
the properties of the simulations that we study. In Section 
|4j first, we present our results in terms of general quantities 
such as the ionization fraction or the optical depth evolu- 
tion. Then, in section [5j we compare the different histories 
of each model enlightened with the merger tree method- 
ology. We finally discuss these results and their limits in 
section [6] and get our conclusions. 



2. Analysing merger trees of HII regions 

2.1. Rationale 

Merger trees have been used for 20 years in e.g. the theo- 
retical studi es of galaxy formation (see e .g. [Lacey fc Cole 
( |1993D and iKauffmann fc White] (|l993| )). More recently 
seminal works such as Furlanetto et al.| ('2004a') and sub- 
sequent analytical and semi- analytical investigations have 
used excursion set formalisms to construct histories of HII 
regions during reionization. Even tough it is strongly con- 
nected to the aforementioned works, the current paper aims 
at using merger trees differently. The trees are extracted 
from simulations^ are used and analyzed as objects worthy 
of studies for their own properties: in other words we an- 
alyze the simulated HII region merger histories from the 
merger tree they produce. The latter is not an intermedi- 
ate to provide higher levels prediction such as e.g. statistics 
on 21cm signal, but is itself a probe of the properties of 
the simulation. On a broader perspective, the methodol- 
ogy presented here could be compared to genus or skeleton 
calculations as a reproductive and quantitative mean to 
discuss the simulated process. It is thus complementary to 
power spectra or time-evolution of average quantities for 
instance. In practice, in addition to tracking individual HII 
regions, we are thus led to analyze the "graph" properties 
of tree-based data as shown hereafter and draw conclusions 
on the simulated physics. Hereafter, we present briefly our 
implementation of the classic FOF technique to identify HII 
regions and the procedure to build merger trees. Additional 
details of these classic tools can be found in the appendix. 



2.2. FOF identification of HII regions 

We have developed a FOF algorithm which is able to sepa- 
rate the different HII regions at a given instant in snapshots 
of simulations and allocate them an identification number. 
The following description of this Friends- of- friends method 
is similar to that used in previous stud ies ( Iliev et al.|2006b 



Shin et al.|2008| [Friedrich et aLpOlT] ). The ionization state 
of the gas is sampled on a regular grid. 



sampled 

First we define the status of the box cells in terms of 
ionization fraction. The rule we have adopted is that a cell 
is ionized if it has an ionization fraction x > 0.5 (as often 
used in other works e .g. [Iliev et al.||2006b[ [Shin et al.|[2QQ8l 



Friedrich et al. 2011). An individual HII region is defined 
by an identification number (ID hereafter). Then if a cell is 
ionized and has an ID, then all its ionized nearest neighbors 
belong to the same ionized region and thus share the same 
ID. Our definition of nearest neighbors is limited to 6 cells 
along the main directions. One single isolated cell is then 
considered as a single ionized region. More details about 
the FOF implementation can be found in appendix [A) 

We have to note that the FOF procedure is not with- 
out caveats. Thereby we could have adopted greater val- 
ues for the ionization fraction threshold for example. We 
can then imagine that HII regions would be more divided 
and more numerous with smaller sizes in this case, possi- 
bly leading to different reionization histories. In appendix, 
we show that the statistics of the regions depend weakly 
on this specific value of the ionization threshold in the 0.3- 
0.9 range, except for the smallest regions with recombining 
episodes. The FOF technique can also connect distant re- 
gions with ionized bridges, an effect that might be undesir- 
able. Alternately we could have chosen other methods for 
the identif ication such as the spherical average method de- 
scribed by'Zahn et al. (2007). In this case the dist ribution 
of HII region sizes would have been smoother (see Friedrich 
[et al.||201l[ ). We have chosen the FOF method because it 
allows to follow regions from snapshot to snapshot in order 
to build the merger tree, which is a feature that the spher- 
ical average method does not allow. We thus have to keep 
in mind the features of the method when we will give our 
conclusion in the next. 



2.3. Merger tree 

The identification of the ionized regions in each snapshot 
of a simulation being known, the merger tree of HII regions 
can be built. Such a tree aims at tracking the identification 
number of an ionized region in order to follow the evolution 
of its properties along time (see appendix [B] for details in 
the merger tree implementation). 

In Figure T] an example of a merger tree is represented 
where each black line represents an ID evolution for a dis- 
tinct HII region. As reionization progresses, there is a grow- 
ing number of ionized regions to follow, as represented by 
the growing number of branches in the tree. We can also 
follow the merger process between HII regions when sev- 
eral branches of the tree lead to the same ID. We observe 
a decrease in the number of regions until there remains a 
single one at the end of reionization. 

From the merger tree thus constructed, we will be able 
to follow the properties of ionized regions during the entire 
simulated reionization. Indeed, it will be possible to study 
the temporal evolution of each HII region, their individual 
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Fig. 1. Illustration of the merger tree of HII regions. Each 
black line represents an ID evolution with the redshift for a 
distinct HII region. For clarity, we represent here only the 
ID evolution for 30 regions. Mergers occur when several 
branches of the tree lead to the same ID. 

merging history, their geometric properties and enlighten 
in a different way the scenario of reionization. 

Let us note that this technique can be sensitive to the 
snapshot sampling of the simulations. Non-regularly spaced 
snapshots or a sparse sampling can weaken the conclusions 
obtained with a merger tree. Preliminary experiments, by 
taking e.g. half of the snapshots, have shown that our re- 
sults seem to be robust, however the quantitative conclu- 
sions may be prone to variations. 

In figure [2] we illustrate the typical properties that can 
be investigated through the merger tree : 

• the number of new HII regions between two snapshots. 

• the number of growing ionized regions. 

• the number of HII regions which recombine. 

• the number of HII regions resulting from mergers. 

• the number of parents involved for a HII region resulting 
from mergers 



3. Simulations 

3.1. Gas Dynamics and Post- processed Radiative transfer 

In this section we aim at describing the cosmological sim- 
ulations of reionization performed and analyzed with the 
merger tree methodology. These simulations were part of a 
se t of experiments at several resolutions and fully described 
in |Aubert & Teyssier| (|2010) and only a brief summary of 
the methodology will be given here. 

The evolution of the gas and sources distribution were 
provided by o utputs of the RAMSES cosmological code 
('Teyssier] |2QQ2 ) that handles the co-evolution of dark mat- 
ter, gas and star particles using an adaptive mesh refine- 
ment strategy. The hydrodynamical equations are solved 
thanks to a second order Godunov Scheme with an HLLC 
Riemann Solver. The gas is assumed to be perfect with 
a 5/3 polytropic index. Metals are included and taken in 
account in the cooling of the gas. The star formation is 
included following the prescription of [Rasera fc Teyssier 
(2006) as well as supernovae feedback. 




Fig. 2. Illustration of the properties that we can follow with 
the merger tree. In each diagram, red items symbolise the 
kind of properties that we follow. 1: the number of new HII 
regions between two snapshots, 2: the number of growing 
ionized regions , 3: the number of HII regions which recom- 
bine, 4: the number of HII regions resulting from mergers 
and 5: the number of parents involved for an HII region 
resulting from mergers 



Initial conditions were gener ated on 1024^ grids, ac cord- 



ing to the WMAP 5 cosmology (Komatsu et al.|(|2009[)) us- 
ing the MPgrafic package and white no ise statistics from 
the Horizon collaboration (Prunet et al. ( 2008) ). Two box 
sizes were considered, 50 and 200 Mpc/h, both with a 1024^ 
coarse grid resolution+ 3 levels of refinement and simu- 
lations were conducted down to 5.5. The refinement 
strategy is quasi-lagrangian with finer levels being triggered 
when the mass within a cell is 8 times the mass resolution. 
The large box is less subject to finite- volume variance ef- 
fects but lacks resolution whereas the 50 Mpc/h box better 
resolve small scales physics (such as star formation) but 
is more prone to cosmic variance and its percolation pro- 
cess during reionization may be already affected by periodic 
boundary conditions 

The radiative transfer was included as a post-processing 
step using the ATON code ( jAubert fc Teyssier 2008 ) in 
all simulations. It tracks the propagation of radiation on 
a 1024^ coarse grid using a moment-based description of 
the radiative transfer equation. It can reconcile the high 
resolution and the intensive time-stepping of the calcula- 
tion thanks to Graphics Processing Units (GPUs) accelera- 
tion obtained through NV idia's CUDA extension to C (see 
Aubert & Teyssier (2010) for further details of the imple- 
mentation). ATON also tracks atomic hydrogen processes 
such as heating/cooling and photo- ionization and can han- 
dle multiple group of frequencies, even though in the cur- 
rent work only a single group of photon has been dealt 
with. The typical photon energy is 20.26 eV, correspond- 
ing to the average energy of an hydrogen ionizing photon 
emitted by a 50 000 K black body: such spectrum is a good 
approximation of a salpeter IMF integrated over the stellar 
particles mass assumed here (see e.g. Back et al. (2009|). 
Calculations were run on 64-128 M2068 NVidia GPUs on 
the hybrid sections of the Titane and Curie supercomputers 
hosted by the CCRT/CEA facilities. Having post-processed 
the simulations, we applied the merger tree methodology 
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Fig. 3. Redshift evolution of the cumulative number of ion- 
izing photons emitted by the sources relative to the number 
of hydrogen atoms. The horizontal line stand for one pho- 
ton per atom. Solid (resp. dashed) lines stand for the 200 
(resp. 50) Mpc/h simulations. Black, blue and red curve 
stand respectively for the 'Star', 'Boosted Star' and 'Halo' 
models of ionizing sources. 



described in the previous section to the ionisation fraction 
field degraded at a 512^ resolution, making them easier to 
analyse. 

3.2. Ionising source models 

At each resolution, three different models were considered 
for the sources: two considered the self-consistent stellar 
particles spawned by the cosmological simulation code. 
The third one uses dark matter halos as proxies for the 
sources. All the models were tuned by trial and error to 
provide a complete reionization by z ~ 6.5 — 5.5 and a 
half-reionization at z ~ 7. Sources were included in the ra- 
di ative transfer calculation f ollowing the same procedure as 
in 



Aubert fc Teyssier (2010). 



3.2.1. 'Star' models 



The st ar formation recipe is described in |Rasera fc Teyssier] 
(|2006 ) and assumes that above a given baryon over-density 
(d ^ 5 in our case), gas transforms into constant mass stars 
(1 X 10^ Moand 2 x 10"^ M© in 200/50 Mpc/h boxes) with 
a given efficiency (e = 0.01). 

The number of stellar particles found at z ~ 8.5 (that 
corresponds to the peak of HII regions number as seen here- 
after) is 8500 (resp. 35 500) in the 200 Mpc/h (resp. 50 
Mpc/h) simulations. Even though it is now standard, the 
modelisation of the formation of 'stellar' particles in cos- 
mological simulations remains a complex and subtle mat- 
ter. Among other effects, it depends strongly on the res- 
olution: the growth of non-linearities is scale-dependent 
and as a consequence the simulated star formation rate 
depends on the ability of the simulation to resolve high- 
density pea ks (see e.g. Springel & Hernquist"(2003),"Rasera" 
Teyssier (^2006) and references therein). Therefore poorly 
resolved simulations tend to develop a population of stars 
at later times and at a slower rate than more resolved ones. 



This limitation is further emphasized during the reioniza- 
tion epoch at large z. As a direct consequence the amount 
of ionizing photons is usually under-evaluated if taken di- 
rectly from star particles even though the situation im- 
proves as the non-linearities evolves in the simulations. For 
instance, with the number of resolution elements and re- 
finement strategy used in the current work, a 12.5 Mpc/h 
box would be necessary to a chieve a convergenc e in th e 
number of emitted photons ( | Aubert & Teyssier ] ( [2qTq| ). 
Otherwise, the star particles may even be too scarce to 
reionize the cosmological box. Furthermore, these sources 
are subject to stochasticity: they contribute to the UV flux 
only during the lifetime of strong UV-emitting stars (that 
we chose to be 20 Myrs). As this component fade away in 
a given 'stellar particle', the latter may not be replaced by 
a neighbour. Until a stationary regime of source renewal 
is installed, these numerical artefacts due to source discret- 
ness and lack of convergence may lead to 'blinking' emitters 
of numerical nature, and therefore to artificial recombining 
HII regions unsustained by inner sources, without any re- 
lation to actual starburst episodes. 

On the other hand, these 'stellar' sources are di- 
rectly generated from the simulation, are simulated self- 
consistently with the gas physics and do not require an 
additional processing to compute the sources locations and 
emissivities. One may suggest that a simple correction may 
be enough to correct from the low resolution effects. The 
simplest is to consider that stellar particles sample correctly 
the location of the main UV sources and a constant fudge 
factor would provide the adequate number of ionizing pho- 
tons to reionize at 2; ^ 7 — 6. Such a model will be named as 
'SXX' hereafter, where 'XX' stands for 200 or 50 depend- 
ing on the box size. A more sophisticated one requires that 
not only the adequate number of photons is generated at 
z = 6 but also at each instant. Since the lack of resolution 
is more critical at early tim es, it implies a stronger c orrec- 
tion at high z than later on. Aubert & Teyssier (2010) have 
shown that the correction is exponential with characteris- 
tic times depending on the resolution. A drawback of this 
type of correction is the fact that sources end up as be- 
ing individually 'decaying' with time: if the renewal rate of 
stellar particles is not high enough, this decay could impact 
the simulated reionization. In this model, 'stellar particles' 
have boosted emissivities per baryon at early times and will 
be named 'SBXX' hereafter, XX standing for the box size. 
It should be noted that the density contrast chosen here to 
trigger star formation can be seen as low and is chosen to al- 
low star formation at our moderate spatial resolution and 
at redshift s where non-linearities are weak at our scales. 
The value chosen here is similar to the one used by e.g. 



Nagamine et al. (2000) at the same resolution and is slightly 
m ore permissive th an most of the thresholds reviewed e.g. 
in Kay et al. (2002) with typical values of ~ 10. However, 
even with such criteria, significant factors should be ap- 
plied to the luminosity of these stellar particles to produce 
enough photons to obtain 'standard reionizations' histories. 
As shown hereafter, the ionized 'patchiness' of these models 
tend to be less structured on small scales than halo-based 
models while nevertheless having their overlap process shar- 
ing the same global properties. Overall it indicates that star 
particles are not overproduced. 

In Fig. [3l one can see that the S200 and S50 'stellar' 
models show different slopes for the cumulative number of 
photons. It results from the fact that the same number of 
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photons must be produced from a smaller number of star 
particles and on a shorter period for the large box. At both 
resolutions, 2 ionizing photons per hydrogen atom have 
been produced when reionization is completed at z ~ 6.5. 
The SB models are as expected 'converged': as the photon 
production become stationary, the same number of pho- 
tons is produced in both box sizes. Compared to the 'Star 
models' emissivities, the 'Boosted' ones are greater at early 
times then smaller for z < 7.5, to achieve the reionization 
at the same epoch. On a more general note, one can see 
that 'stellar' particles appear earlier at high resolution as 
expected. 

Quantita t ively, we followed the simple modelisation of 



Back et al. (2009) to assign an emmissivity of 90 000 
photons per stellar baryons over the lifetime of a source. 
For S50 and S200 the enhancement factors to account for 
convergence issue are respectively 3.8 and 30, while the 
boost temporal evolution of the SB50 and SB200 is equal 
to max(l,aexp(A:/t)) with (a, A:[Myr]) = (1.2,1500) for SB50 
and (1.2,3000.) for SB200, t being the cosmic time. 
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Fig. 4. Evolution of the volume weighted average ionization 
fraction with redshift for the six simulations. 



3.2.2. Halo sources 



As a second choice we adopt a semi-analytical model for 
the production of ionizing sources based on dark matter 
halos. Each halo is assumed as a star formation site with 
an emissivity proportional to th e halo mass . This p rocedure 
is largely inspired by the work of Iliev et al. (2006b) where a 
constant mass-to-light ratio is assumed such as the ionizing 
flux of each halo has the following expression: 



= aM 



(1) 



where a is the emissivity coefficient. Here we have 
chosen values of a = 5.9 x 10^^ and a = 3.5 x 10^^ 
photons/s/M© for both 200 and 50 Mpc/h boxes respec- 
tively. These values were chosen in order to reach a reion- 
ization at z ~ 6.5 — 6 similar to the 'Star' and 'Boosted 
Star' mo dels. Halos were detec ted using the parallel FoF 
finder of Courtin et al. (2011), with a linking length of 
b=0.2 and a minimal mass of 10 particles. This minimal 
mass corresponds to 9.8 x lO^M© for the 50 Mpc/h box 
and 6.3 x lO^M© for the 200 Mpc/h box, implying that all 
halos can be considered as emitters. The number of halos 
found at z - 8.5 is 33 100 (resp. 285 000) in the 200 Mpc/h 
(resp. 50 Mpc/h) simulations. It should also be emphasized 
that, contrary to 'Star' and 'Boosted Star' models, a whole 
spectrum of mass, and thus emissivities, is available for ha- 
los. No lifetime has been assigned to these sources and a 
halo would produce photons as soon as it is being detected 
until it disappears or merge. In the following sections we 
will refer to this semi-analytical star formation simulations 
with the following acronyms: "H200" or "H50" for both box 
sizes of 200 and 50 Mpc/h. 

In Fig. [3j one can see that 'Halo 50' follows closely 
'Star 50', implying that provided the appropriate correc- 
tion, stars can mimic halos emission history. A greater dif- 
ference can bee seen in the 200 Mpc/h box, where the 'Halo 
200' model shows a greater slope than the 'Star 200' model: 
there exists a greater difference in the buildup of halos com- 
pared to the star population at lower resolution, due mostly 
to the difficulty to achieve star formation density thresholds 
in large scale simulations. At z ~ 6.5, roughly two photons 
per hydrogen atom were produced in all cases. 



Boosted star 200 
Halo 200 
Star 200 
Boosted Star 50 
Halo 50 

Boosted Star 50 
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Fig. 5. Evolution of the Thomson optical depth with red- 
shift for the six simulations. 



4. Global features 



4.1. lonisation fraction 



Figure [4] shows the evolution of the volume weighted aver- 
age ionization fraction as a function of redshift for the three 
models and for both 200 and 50 Mpc/h boxes. Average ion- 
ization curves show a similar value of<x>'^0.5ata red- 
shift of z ^ 7. At first glance we see that every models are 
comparable in terms of ionization fraction. The boosted 
star model just shows an earlier "take off' in the ioniza- 
tion curve than in both others models for the two box 
sizes. Conversely in the Star and Halo models the ioniza- 
tion curves "take off' later {z ^ 10 — 9) but the reionization 
is achieved earlier than in the boosted Star model. These 
differences are attenuated in the 50 Mpc/h box where the 
ionization histories become very close from one model to 
another. 
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4.2. Optical depth 

Figure [5] shows the evolution of the optical depth as a func- 
tion of redshift for the six simulations given by 



t{z) cat 



J z 



At 



(2) 



where Gt is the Thompson cross section of the electron and 
ne{z) =< x{^z) > nniz) is the density of electrons released 
by ionized hydrogen atoms at redshift z. We also represent 
the constraint range obtained from the five year release 
of CMB measur e ment made by the WMAP collaboration 
dKomatsu et al.| ( |2QQ9| )) at the 1 a level. 

We immediately see that all the simulations converge 
in terms of optical depth. The 200 Mpc/h box simulations 
reach a same value at z ~ 8. The only difference being that 
the optical depth is slightly greater in the boosted Star 
model before z ^ 8. This is naturally explained by the fact 
that the ionization history is more extended in this model as 
seen with the ionization curve in figure |4] In the 50 Mpc/h 
box, the curves become almost superimposed from z ^ 10. 
Again this comforts us in saying that all the simulations 
are comparable. 



4.3. lonisation fields 

The maps of figures [6] and [7| show the ionization field maps 
for the three models for both 200 and 50 Mpc/h boxes. 
In each figure, colour encodes an individual HII region de- 
tected by the FOF procedure. 

At z ^ 10.5 and z ~ 8.4, the basic features of the 3 
type of reionization can be spotted: SB models exhibit a 
few regions with large radii, whereas H models have smaller 
and numerous individual regions. S models are intermediate 
with larger ionized patch than detected in the halo-based 
simulations but more individual regions than found in the 
SB models. As such it refiects the differences in source mod- 
eling where halos are more numerous than stellar particles 
and therefore share ionizing photons over a larger number 
of weaker sources. On the other hand, SB models produce 
large regions at the earliest redshifts because of the strong 
initial correction to the source emissivity and which favors 
large regions and potentially early overlaps. 

At the end of the reionization at z ~ 6.9 and in all mod- 
els, a single large HII region is detected by the FOF method 
that results from a connected network of multiple ionized 
regions. Nonetheless, the Halo model presents a greater re- 
silience to percolation as it presents a map more structured 
with more individual HII regions. 

When the spatial resolution is increased with the 50 
Mpc/h box, the figure j?] shows the same tendencies as in 
the 200 Mpc/h box. The H50 map presents the most im- 
portant number of regions with the smallest sizes until late 
phases of the reionization and the SB50 still presents the 
lowest number of regions with the largest sizes. S50 still 
represents an intermediate case. Overall and as expected by 
the greater spatial resolution, ability to capture the clus- 
tering and larger number of sources in a 50 Mpc/h box, 
these maps present a greater level of granularity and more 
individual regions than the 200 Mpc/h versions. 



5. Merger tree properties 

Our main interest is to investigate the impact of the ionizing 
source models on the reionization history. Typically, we will 
see what quantities are retained from one model to another 
and how these models imprint differences on the history 
observed. 



5.1. Evolution of the number density of HII regions 

The figure [S] shows the evolution of the number density of 
ionized regions as a function of redshift for the three models 
of ionizing sources and for both box sizes. We distinguish 
in addition the number of each type of HII regions: new re- 
gions, those resulting from mergers, those expanding with- 
out merging and those recombining are respectively shown 
in blue, green, red and yellow. We note that the distri- 
butions present the same general shape regardless of the 
source model. They all show their maximum number of re- 
gions at broadly the same redshift z peak for both spatial res- 
olution respectively at z ~ 8 and z ~ 9 for 200 Mpc/h and 
50 Mpc/h. Before Zpeak, HII regions appear and expand, 
populating the box with more and more regions. We can 
refer to this period as the " pre-overlap" period where the 
reionization is dominated by the expansion and the birth 
of new HII regions. Conversely after Zpeak, HII regions be- 
gin to merge intensively and decrease in number during 
an "overlap period", that continues until the reionization 
is completed. Even though the shape of the distributions 
are similar, the absolute number of HII regions is much 
greater in the Halo model than in the two 'stellar ones' 
and with a greater discrepancy in the 50 Mpc/h box. This 
absolute number remains within the same order of magni- 
tude between the Star and the boosted Star model for both 
boxes. This difference is not surprising and typically refiects 
the difference in the relative number of ionizing sources 
produced from one model to another. Indeed the number 
of dark matter halos assumed as ionizing sources in the 
semi analytical model is greater than the number of auto- 
consistent ionizing sources generated in both Star models : 
at z=8.5 there are ~ 4 (resp. 8) times more halos than 
stellar particles in the 200 Mpc/h (resp. 50 Mpc/h) and it 
corresponds to the factor measured in terms of detected HII 
regions. Finally, when splitting these distribution in differ- 
ent types of ionized regions (new, expanding, merging and 
recombining), they follow the overall pattern of increase 
and decline after z peak with differences detailed in the next 
section |5.2[ Let us only mention the case of new regions 
that also follow the same pattern even though at face value 
mergers only affect pre-existing ionized patches. It should 
be reminded that as reionization progresses the amount of 
neutral volume decreases, limiting the possibility of having 
new sources in a neutral area. Hence, new regions are also 
affected by the overlapping process, albeit indirectly. 



5.2. Number and volume fraction of the different types of HII 
regions 

To emphasize differences in the histories of the various kind 
of regions. Fig. [9] shows the evolution of their relative pro- 
portion (rather than their number density) in terms of num- 
ber and Fig. [lO] presents the volume fraction of the whole 
simulated boxes occupied by each type of HII region. 
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Fig. 6. lonisation map for the three ionizing sources models for three distinct redshift for the 200 Mpc/h size box. The 
colours encode the different identification number allocated to the HII regions by the FOF algorithm 



5.2.1. Recombining regions 

The most striking feature of SB models (and to a lesser 
extent S models) is the presence of recombining HII re- 
gions, tracked as connection-less branches in the merger 
tree. Recombination can be driven by two main factors: 
high density features or evolving sources. The evolution of 
the density is common to all experiments and the lack of 
such regions in H models indicates that it is unlikely to 
be the main origin of these peculiar regions. Nevertheless, 
a hint for this 'density' effect can be seen in the larger 
amount of such regions in the S50 model compared to S200, 
or their marginal presence in H50: higher resolution leads 
to higher contrasts and thus higher recombination rates. 
On the other hand, stellar sources are more prone to vari- 
ation than 'halo' sources. First they exhibit some level of 
stochasticity, hence sources that turn off may not be re- 
placed through an efficient renewal process in order to sus- 
tain some regions, whereas halo sources have by construc- 
tion a continuous emission. It explains the lack of such re- 
gions in H models and their detection in S and SB mod- 
els. Furthermore SB models have sources with a decreasing 
individual emissivity, leading to a typical scenario where 
strong early sources produce large regions that see their in- 



ner ionizing engine becoming progressively weaker, or stars 
being replaced by weaker ones, eventually leading to re- 
combinations. Combined to stochasticity, it can lead to re- 
gions that are being 'dissolved', potentially at several places 
within one region. Increasing the resolution tend to dimin- 
ish the contribution of this type of regions in SB50 models: 
the differential boost applied to source is weaker than in 
SB200, hence decreasing the effect of inner exhausted en- 
gine and the stellar renewal is more efficient, leading to a 
production of ionized regions closer to a stationary regime. 
Looking at Fig. [lO} it is interesting to note that in all the 
cases where they are detected, recombining regions stand 
for a small fraction of the total ionized volume: between 
^ 10~^ and ~ 10~^ of the total volume of the box for the 
H50 and the SB200 model respectively. Therefore, even if 
recombining regions are good indicators of the behaviour of 
the source models, it can be seen that they correspond only 
to an almost negligible fraction in terms of volume and can 
be considered to some extent as marginal. 
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Fig. 7. lonisation map for the three ionizing sources models for three distinct redshift for the 50 Mpc/h size box. The 
colours encode the different identification number allocated to the HII regions by the FOF algorithm 



5.2.2. Regions in expansion 



Expanding regions are detected as having one 'parent' only 
and distinguish HII regions that expand without additional 
event. In absolute number (see Fig. |8| their number peaks 
at 2;peak and their proportion is always dominant until the 
very end of the reionization. Incidentally, it shows that the 
temporal sampling is high enough to generally track the 
initial expansion of individual regions. Also, the fact that 
such regions are still detected at the later stages suggests 
that even at the later stages of the percolation, not all the 
regions are involved in a global merging process: some re- 
gions still have room to experience stages of 'quiet growth' 
before eventually be part of the global ionized background. 
In terms of volume (see Fig. 10), expanding regions are 
the main contributors to the total ionized volume in every 
models until typically ^peak- This is not surprising since 
the expanding regions are naturally much greater in size 
than the new regions and are always dominant in number 
(see Fig. [s]). Once the overlap start to be efficient, merged 
HII regions superseed the expanding ones in volume, as ex- 
pected. 



5.2.3. New regions 

New regions (without parents in the tree) track the for- 
mation of new sites of emission, i.e. star or halos forming 
in areas distant enough of pre-existing ionized volumes. Of 
course, any source that appears within an ionized region 
because of clustering won't contribute to this population. 
First, their proportion (and their absolute number) present 
regular 'spikes'. The tree is constructed using a time sam- 
pling different than the sampling used to include sources 
in the radiative transfer calculation: the tree sampling fre- 
quency is low enough that between two snapshots at least 
a new generation of source has been included and new 
regions are thus detected. But this frequency is also low 
enough that regularly, two generation of new sources were 
included in the RT calculation between two snapshots, cre- 
ating bursts of new regions. The presence of these spikes 
does not affect our analysis and conclusions. 

Globally all experiments show an initial gradual de- 
crease of the proportion of new regions. In the 50 Mpc/h 
experiments, this decline stays until the end of the reioniza- 
tion. It is related to a pile-up effect of pre-existing HII re- 
gions: until 2;peak their number increase at a higher rate than 
the one of new regions. For instance, regions keep growing 
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Fig. 8. Evolution of the number density of each kind of HII regions as a function of redshift for the three models of 
ionizing sources and for both boxes of 200 and 50 Mpc/h. Panels (a), (b) and (c) respectively represents the distribution 
for the Boosted Star, the Star and the Halo model for the 200 Mpc/h box, while the panels (d), (e) and (f) are for the 
same models but for the 50 Mpc/h box. The colours stand as follow: Blue for the new HII regions. Red for the expanding 
regions. Yellow for the regions that will recombine and Green for the regions resulting from mergers. 



without efficient merging over several generation of new HII 
regions: it increases the number of pre-existing regions and 
decrease mechanically the proportion of the new ones at a 
given moment. After 2;peak, the proportion of new regions 
increases in 200 Mpc/h experiments and arguably remains 
constant in 50 Mpc/h ones. By construction the overall 
number of ionized sites decreases in this 'post-overlap' pe- 
riod, indicating that among these sites, new regions suffer 
to a lesser extent of the percolation process. Hence while 
the destruction rate of pre-existing HII regions is impor- 
tant, new stars or halos manage to appear in neutral areas, 
sustaining (in the 50 Mpc/h box) or even increasing (in the 
200 Mpc/h box) the relative contribution of the new re- 
gions. This increase of the contribution of new regions seen 
in S200 and H200 is indicative that the destruction of pre- 
existing regions is more sudden than at higher resolution: 
it was hinted in Fig. M where the post ^eak evolution is 
smoother at high resolution or in Fig. |4] where the average 
ionized fraction exhibits a sharper evolution at reionization. 
Regarding this sharper evolution, let us emphasize that a 
more accurate tuning of the source emissivity could have 
increased the matching of the reionization history of the H 
and S models and may be the origin of this difference be- 
tween high and low resolution. It may also be the result of a 
different behaviour in the propagation of fronts at low res- 
olution (with faster fronts and weaker shielding), leading 
to a more radical percolation in the 200 Mpc/h simula- 



tions than in the 50 Mpc/h boxes. Finally, let us mention 
that the presence of a strong recombining component in the 
SB200 model lead to a totally different evolution of the new 
regions contribution: it remains broadly constant at all red- 
shifts with a weak dip at z ~ 10. Without this recombining 
regions, it would also exhibit a rising relative weight of new 
regions, albeit starting much earlier than %eak- 
Regarding their volume fraction (see Fig. 



10), these 



newly detected regions only represent a maximum of ^ 
10~^ — 10~^ of the total volume even if there is a siz- 
able number of this kind of regions detected. This is un- 
derstandable because of their very nature: new regions are 
by definition expected to be small, since they just ap- 
peared and their number cannot compensate this effect. 
Interestingly, the total volume occupied by the new regions 
remains broadly constant until overlap : it indicates that 
even though sources increase in number with time (and thus 
the potential number of new regions) it competes with the 
gradual decrease of available neutral volume, these two ef- 
fects 'conspiring' to keep constant the volume imprint of 
the new patches on the network of HII regions. 



5.2.4. Regions resulting from mergers 

These regions have more than one parent and result there- 
fore from the merger of several regions. They track the co- 
alescence of pre-existing regions into larger ones and even- 
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Fig. 9. Evolution of the proportion of each kind of HII regions as a function of redshift for the three models of ionizing 
sources and for both boxes of 200 and 50 Mpc/h. Panels (a), (b) and (c) respectively represents the distribution for the 
Boosted Star, the Star and the Halo model for the 200 Mpc/h box, while the panels (d), (e) and (f) are for the same 
models but for the 50 Mpc/h box. The colours stands as follow: Blue for the new HII regions. Red for the expanding 
regions. Yellow for the regions that will recombine and Green for the regions resulting from mergers. The black vertical 
line shows the peak of the absolute number of HII regions: Zpeak- 



tually track the final overlap into a single large HII region 
at the end of reionization. This last stage can bee seen in 
all models where there exist a late snapshot where 100% 
of the detected region result from mergers, an indication of 
the ultimate merger. S and H model present similar evolu- 
tions for the merger populations. Their proportion is 'co- 
incidentally' maximal at Zpeak- as the number of ionized 
patches gets larger, a larger fraction of them is involved in 
mergers, indicating a 'crowding' or clustering effect where 
a smaller volume is available for expansion as the number 
of HII regions increase. Later on, for z < Zpeak, the propor- 
tion of mergers decreases indicating that only a subset of 
the ionized regions do actually merge, potentially only one 
that would phagocyte the others, reducing the overall num- 
ber of individual regions until the end. Hence the peak of 
mergers fraction could be seen as the rise of one or several 
dominant regions and this rise appears when the number 
of individual region is maximal. 

Comparing the two resolutions, it can be noted that 
mergers are concentrated over a smaller range of redshifts 
at low resolution whereas merger can be detected at a sig- 
nificant level during the whole experiment in the S50, H50 
and SB50 models. Firstly, it can be the consequence of the 
slightly more extended history of reionization in 50 Mpc/h 
simulations, already mentioned regarding the {x{z)) trends 
(see Fig. |4|. Furthermore, sources and individual regions 



are more numerous at any time and in a smaller volume, in 
favor of a more generalized contribution of mergers over a 
large range of redshifts. 

Finally, it should be noted that the SB200 model 
presents an early peak in the fraction of merger, at z ~ 10 
instead of z ^ S for all the others experiment, with an al- 
most zero contribution later on, even though the absolute 
number of regions manage to decrease at some point. It 
suggests that the onset of the global percolation started 
earlier, due to early large HII regions induced by the large 
initial boost of this model, and created a single region that 
monopolize the merger process. Interestingly, in this case 
Zpeak occurs later, indicating that while this main region 
grows, there remains a significant amount of neutral vol- 
ume to host the apparition of new regions. 

Finally, Fig. [lO] clearly indicates that these merger re- 
gions dominate the fraction of the ionized volume for z < 
Zpeak- It was uot obvious given that their number density is 
always the lowest contribution (cf. Fig. [8|. It is indicative 
that few of these regions produce a network that dominate 
the ionized volume and/or that such regions are individu- 
ally very large, as one would expect from them being the 
result of several mergers. This point is assessed in greater 
detail in the next sections. 
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Fig. 10. Evolution of the volume fraction of each kind of HII regions as a function of redshift for the three models of 
ionizing sources and for both boxes of 200 and 50 Mpc/h. Panels (a), (b) and (c) respectively represents the distribution 
for the Boosted Star, the Star and the Halo model for the 200 Mpc/h box, while the panels (d), (e) and (f) are for the 
same models but for the 50 Mpc/h box. The colours stands as follow: Blue for the new HII regions. Red for the expanding 
regions. Yellow for the regions that will recombine and Green for the regions resulting from mergers. The black vertical 
line shows the peak of the absolute number of HII regions: Zpeak and the black horizontal line shows the volume fraction 
of one cell of the grid. 



5.3. Sizes of HII regions 

5.3.1. Sizes distribution with redshift 



Fig. 11 shows the evolution of the distribution of the HII 



regions radii at each instant as a function of redshift for 
the 200 and 50 Mpc/h box. Like Friedrich et al.| ( |201l[ ), we 
compute the volume of the HII region and then we derive 
the effective radius corresponding to a sphere of equal vol- 
ume through the following expression R = [3/(47r)V]^/^. 
The white dashed curve shows the radius evolution of the 
last HII region which remains at the end of the simulation. 
With the help of the merger tree we follow this region back 
in time and calculate the radius of its main progenitor at 
each instant. In addition the solid white line represents the 
evolution of the average radius of HII regions as a function 
of redshift. The color code in the distributions is made so 
that the brightest red cell contains only a single HII re- 
gion. Other red tones up to the brightest green designate 
cells populated by a couple HII regions up to a value cor- 
responding to half the maximal value of the distribution. 
Finally the purple tones (from the darkest to the brightest) 
encode the largest values of the distribution. 

Typically, each distribution seems to trace the under- 
lying ionizing source prescriptions related to the model 
considered. The Star model has a constant star emissiv- 



ity with redshift. This is reflected in the radii distributions 
where ionized regions are concentrated with a constant ra- 
dius about r ~ 4 Mpc/h over the whole range of redshift in 
the 200 Mpc/h box. We find a gap at about r ~ 2 Mpc/h in 
the size distribution where there are only a few regions with 
radii under this constant radius. This would indicate that 
the time sampling used here cannot entirely capture the 
fast tracks followed by the HII regions in the radius-redshift 
space with strong inner sources. This effect is attenuated for 
the 50 Mpc/h box (that contains typically weaker sources 
and HII regions with smaller growth rates) where we ob- 
serve a more continuous distribution even if there remains 
a typical cutoff radius of r ~ 0.4 Mpc/h, with smaller oc- 
cupation numbers at low values. 

On the other hand, the Halo model implies that each of 
them has an emissivity proportional to its mass. As these 
masses cover a large range, we find a large range of radii 
for the resulting HII regions in the distribution of the 200 
Mpc/h box. The regions are concentrated in larger inter- 
vals of radii than in the other models that typically trace 
the underlying mass range of halos. The shape of the distri- 
bution is almost the same when we consider the 50 Mpc/h 
experiment, with smaller regions as expected since smaller 
halos are available at higher resolution. 
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Fig. 11. Evolution of the HII regions radius distribution as a function of redshift for the three models of ionizing sources 
formation and for both boxes of 200 and 50 Mpc/h. Panels (a), (b) and (c) respectively represents the distribution for 
the Boosted Star, the Star and the Halo model for the 200 Mpc/h box, while the panels (d), (e) and (f) are for the 
same models but for the 50 Mpc/h box. The color code stands as follows: the brightest red cells represent the location 
in the distribution populated by a single HII region, other red tones up to the brightest green cells spans distribution 
densities from a couple HII regions up to a value corresponding to half the maximal value of the distribution. Finally the 
purple tones (from the darkest to the brightest) encode the maximal values. The evolution of the average radius and the 
evolution of the radius of the main region are shown respectively with the solid and dashed white lines. 



Finally, the boosted Star model has a boost for the star 
emissivity that decreases with redshift. Indeed the distri- 
butions show, at early times, large HII regions detected 
without small counterparts. This gap would be the result 
of the very powerful boost for ionizing sources at high red- 
shift combined to the time sampling of the simulation which 
allows to detect only regions when they have a large radii 
early in the reionization. Then we can observe a decreasing 
gradient for the typical radius of the regions as the boost 
for the emissivity decreases with time. Surprisingly we find 
a bimodal distribution from a redshift of 2; ^ 9 until the 
reionization at 2: ~ 6. Some HII regions are concentrated 
with radii below or under ~ 1 Mpc/h with a gap in the 
distribution for this value. To a lesser extent, the same ef- 
fect could be seen in the Star 200 model. This regions with 
radii under ~ 1 Mpc/h could be the recombining regions 
that we have found in figure |9] This could also be com- 
bined with the fact that around these redshifts the boost 
emissivity falls down under a certain value. Then, the time 
sampling of the simulation would become fine enough to 
detect some new or expanding regions with these radii. We 



will further investigate this in the following section 5.3.2 



For the SB50 model, this bimodal distribution disappears: 
a weaker boost amplitude and variation combined to a more 



stationary production of stars promotes a better tracking 
of regions sizes as they grow and reduce the contribution of 
recombination, as seen earlier. 



Considering again the S and H models, for both box 
sizes, it is interesting to note that the S model produce a 
truncated version of the H one. Above the cutoff radius (0.4 
Mpc/h (resp. 4 Mpc/h) in the S50 (resp. S200)) the S and 
H distributions are quite similar. At smaller radii, S has less 
objects than H. It is indicative that there is a scale above 
which the clustering of halos combined to their own mass- 
proportional emissivities produce a similar radii distribu- 
tion to the one provided by the clustering of stellar sources 
and their own constant emissivity. At smaller scales, the 
stellar sources are too scarce to reproduce the regions cre- 
ated by halos, resulting in overpowered individual sources 
with large radii. At larger scales the two approach produce 
equivalent size distributions and with an appropriate cal- 
ibration, hydrodynamically created sources can match the 
halo result. It should be noted that this typical scale ap- 
pears at smaller values in S50 models as sources are more 
numerous and are less prone to stochasticity and therefore 
converges toward the halo model behavior on smaller vol- 
umes. 
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Finally, all the experiments present a single main HII re- 
gion that dominates in size. It appears at 2: ~ 8 in the S200 
and II200 models and much earlier {z ~ 9) in the SB200. 
As suggested by the previous study on merger populations, 
it is likely the consequence of the boost that creates early 
large regions that merge at high z to give the dominant one. 
Interestingly the dominant region of the late reionization is 
not dominant in size at every redshift: this is indeed the 
case for the II200 model and to some extent for the S50 
model too, but in all the other case this region started as a 
non special one inside the population, at least from the 'ra- 
dius' point-of-view. Furthermore, the first progenitor of the 
dominant region is among the very first ionized regions of 
the II200 and II50 models but appears later in the S and SB 
experiments. It is noteworthy that even though we argued 
that S and H model produce similar populations (for regions 
larger than the cutoff radius) it seems that for individual 
cases, differences can persist, especially for the buildup of 
the dominant region. As we will see in next sections, the 
channel through which this specific region grows (through 
expansion rather than merger) also differs as hinted by the 
difference between the smooth evolution of the dominant 
region of H models and the sharp, kinked and late rise of 
its equivalent in e.g. the S200 model. 

5.3.2. Radius distribution of the different kind of HII regions 



Fig. 12 presents the size distribution for each kind of HII 
regions as a function of redshift for all three models and 
for both 200 and 50 Mpc/h boxes. Here we show the same 
distribution seen in figure [TT] but by plotting the contours 
related to each type of HII regions that we can discern. 
The colour code stands as follows: Blue corresponds to the 
new HII regions, red to the expanding ones, green to those 
resulting from mergers and yellow to regions that will re- 
combine. Once again we show the evolution of the average 
radius of the region with the black curve and the evolution 
of the radius of the last HII region with the black dashed 



curve as in the figure 11 



Initially, we observe that each kind of HII region oc- 
cupies a dedicated range of radii in the whole distribution 
for all models and for both 200 and 50 Mpc/h box sizes. 
As expected new regions occupy the smallest range of radii 
while the regions resulting of mergers are those which pop- 
ulate the top of the distribution. The expanding HII regions 
are in the middle with radii greater than the new ones and 
smaller than those that will merge. We can also see in each 
model that the peak of the distribution of the figure [TT] cor- 
respond to a radius range that is simultaneously covered by 
the three kinds of HII regions. Alternatively, we could have 
imagined a clear separation between the different range of 
radii covered by the different kind of HII regions instead 
of the observed overlap of the different distributions. This 
tells us first that the time sampling here allows us to detect 
HII regions with same radii but belonging to different kind 
of HII regions with different paths in the radius-redshift 
space. It also indicates that the peak of the radii distribu- 
tion must be understood as the most likely radii of detec- 
tion, with a mix of HII regions at different stages of their 
evolution (newborn, growing or merging). 

Interestingly, for each kind of regions, we find the prop- 
erties of the underlying ionizing source models in the evolu- 
tion of the radius distribution with redshift. In other words 
the R(z) distribution for the size evolution of each kind of 



regions as function of redshift is dictated by the R(z) dis- 
tribution of the new regions and is propagated to the other 
kind of regions. We thus can see a similar R-z relationship 
for the expanding region and those resulting from merger 
with a simple translation. In the Star models S50 and S200 
the new regions appear with a relative constant range of 
radii during the whole reionization; again, it reflects the 
typical constant emissivity of the stars appearing during 
the whole simulation combined to the time sampling used 
here. The expanding regions and those resulting from merg- 
ers also show this constant radius range but translated to 
greater radius ranges. The S50 allows to track smaller new 
regions thanks to its higher resolution in terms of sources 
with smaller individual emissivity: as a consequence all sub- 
sequent types show smaller typical radii than in the S200 
model. 

The H200 model presents new regions that appear with 
sizes more distributed than in the S50/S200 models and 
typically result from the mass range of the dark matter ha- 
los that are assumed as ionizing sources. Then the expand- 
ing regions cover a much greater radius range compared 
to the two other models and we find the same trend for 
the regions resulting from mergers. The last point indicates 
that the merger process is not dominated by a few large 
regions as in the two other models but also involves sets 
of small HII regions. In the H50 model the mixing is more 
pronounced and the ranges of radii covered by each kind 
of regions overlap. Highly clustered regions are more rep- 
resented at high resolution where dense set of small halos 
with small regions can also merge. In H200, mergers occur 
at later stage of radii evolution, in a more homogeneous 
distribution of halos at low resolution. 

Finally, SB200 shows that new regions appear with de- 
creasing radii as reionization progresses and as the boost for 
the star emissivity decreases. The expanding regions and 
those resulting from merger show the same R-z decreasing 
trend, seeded by the new regions population properties. It 
should be noted that at some point {z < 9), new regions 
can be detected at small radii: the declining boost leads to 
individual regions that appear with smaller size and growth 
rate, increasing the likelihood to spot them at small radii. A 
weaker R-z slope is found in the SB50, as a consequence of a 
smaller boost correction at higher resolution. As expected, 
new regions with small radii are found earlier (for z < 12) 
and globally share a greater number of broad features (like 
the range of radii of each type) with the S50 model. 

Let us finally focus on the HII regions that recombine 
in the size distributions of each model. We find as expected 
in section [5.3. 1| that regions that will recombine are those 
with radii below the radius of ~ 1 Mpc/h in both SB and 
S models and at both resolution. As already mentioned, 
the recombination is driven either by density or a modifi- 
cation of the sources. In the case of SB models, it is fairly 
obvious that strong early sources create large HII regions 
that may be unsustainable by later generation of sources 
which are individually weaker by construction of the boost. 
Furthermore, the boost guarantees a global convergence but 
locally a strong early source can be replaced by too few 
weaker later sources, especially if the source renewal is sub- 
ject to stochasticity. Once sources become too dim, it leads 
to the fragmentation of regions and recombinations at small 
scales. S models also present some degree of recombination, 
even though they do not contain declining sources. These 
recombinations can be explained by local dense patches and 
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(d) Boosted Star 50 (e) Star 50 (f) Halo 50 

Fig. 12. Evolution of the radius distribution of each kind of HII regions as a function of redshift for the three models 
of ionizing source formation and for both boxes of 200 and 50 Mpc/h. Panels (a), (b) and (c) respectively represents the 
distribution for the Boosted Star, the Star and the Halo model for the 200 Mpc/h box, while the panels (d), (e) and (f) 
are for the same models but for the 50 Mpc/h box. The colours stand as follows: Blue for the new HII regions. Red for 
the expanding regions. Yellow for the regions that will recombine and Green for the regions resulting from mergers. In 
addition the solid and dashed black lines represent respectively the radius evolution of the main region end the evolution 
of the average radius for the HII regions. 



it could be supported by the fact that such regions only ap- 
pear at later time {z < 9) in the S200 model, when local 
clustering is high enough. Furthermore, in both S50 and 
SB50 models, these regions can be found over extended pe- 
riod compared to their 200 Mpc/h equivalent, related to 
the fact that dense patches that trigger recombination are 
more likely to be found at high resolution. 



5.4. Mergers of HII regions 

5.4.1. General evolution of the merger process of HII regions 



Fig. 13 is restricted to regions that have more than one par- 
ent, i.e. that result from mergers. It presents the redshift 
evolution of the distribution of the number of parents. It 
can be seen as a measure of the patchiness or granularity 
of the overlapping process and is also a first indication of 
the inner complexity of the regions: an ionized volume that 
result from tens of progenitors is shaped (like e.g. the front 
propagation, its growth rate, or the inner distribution of UV 
flux) by the different properties of at least tens of sources 
whereas an HII region coming from a binary merger is likely 
to be, maybe naively, more straightforward to relate to its 
inner sources. In addition, the solid white line in flgure [13] 
shows the evolution of the average number of parents and 



the white dashed line stands for the evolution of the num- 
ber of parents for the last remaining region at the end of 
the simulation. Again, with the help of the merger tree we 
follow this region back in time and evaluate its number of 
parents at each instant. The color code in the distributions 
is identical to the one used in the flgure pT 



Firstly, in most of the models the mergers between re- 
gions operates in a binary-tertiary manner as seen with the 
peak at ~ 2 — 3 in the distributions. It is essentially an indi- 
cation that the temporal sampling of the tree is flne enough 
to track individual mergers whereas a coarser grained sam- 
pling would have presented a typical larger value. One ex- 
ception is the H50 model, where clearly such mergers are 
missing from 2; ~ 10 to z ~ 7. In this case, the abundance of 
clustered individual regions at high density is much higher 
than in other models and leads to multiple mergers during 
a single time step for a given region. 

Then, in every model we flnd the emergence of a sin- 
gle region which is the result of more mergers than other 
regions. In the S200 and H200 models, this region appears 
at z ~ 8 and earlier in the SB200 one, at z ~ 9. The same 
discrepancy can be found at higher resolution with a forma- 
tion at z ~ 9 for the S50 and H50 models and 2: ~ 10 in the 
SB50. Interestingly, this merger-dominant regions appear at 
the same moment as the radius-dominant region discussed 
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Fig. 13. Distribution of the number of parents for the HII regions resulting from merger as a function of redshift for the 
three kinds of ionizing sources and for both box sizes of 200 and 50 Mpc/h. Panels (a), (b) and (c) respectively represents 
the distribution for the Boosted Star, the Star and the Halo model for the 200 Mpc/h box, while the panels (d), (e) and 
(f) are for the same models but for the 50 Mpc/h box. The color code stands as follows: the brightest red cells represent 
the location in the distribution populated by a single HII region, other red tones up to the brightest green cells spans 
distribution densities from a couple HII regions up to a value corresponding to half the maximal value of the distribution. 
Finally the purple tones (from the darkest to the brightest) encode the maximal values. We represent in addition the 
evolution of the average number of parents and the evolution of the number of parents of the main region respectively 
with the solid and dashed white lines. 



in the previous section. In Fig. [Tsj the white dashed line 
shows the evolution of the number of parents of the radius- 
dominant region and unsurprisingly it matches the path of 
the merger-dominant one. Thus, a moment exists when an 
HII region begins to monopolize the merger process in the 
box. Then the region grows faster by successive mergers and 
will quickly dominate other regions in terms of size. From 
this instant (called zbkg)^ the probability that a location 
is part of the global UV background instead of being irra- 
diated by a local source increases and can be considered as 
the onset of the overlapping process. This is typically the 
moment when on average the local information about the 
reionization process start to be lost to the benefit of a single 
reionization history. 

Size is an important factor for regions to get high 
numbers of progenitors as a larger volume naturally pro- 
mote encounters with distant ionized patches. However it 
can be noted that the radius-dominant region is only one 
among the merger-dominant regions. Clearly the distribu- 
tion shows regions with larger number of progenitors and 
this is especially the case for H200 and H50 models: the 
high density of sources and their clustering can induce high 
merger rate in smaller volumes. More generally, high reso- 



lution experiments present greater merger rate thanks to a 
higher number of sources and higher local clustering. Also 
H models can achieve number of progenitors that are ten- 
fold greater than the S and SB equivalent as another con- 
sequence of large number of halos. The two stellar models 
differ also in the number of regions with a large number of 
progenitors: close to zbkg, S models present clear detection 
of regions with a number of progenitors close to those found 
for the radius-dominant one (10-50 progenitors), whereas in 
SB models such regions are hardly found. It indicates that 
the early buildup of a dominant region in the SB model 
tends to prevent the formation of separate regions that ag- 
gregate large number of sources on their side, probably be- 
cause they were incorporated early in this dominant vol- 
ume. At higher resolution the discrepancy is weaker, but 
in the S50 models, the dominant region share with more 
other regions the property of having a large (> 10) number 
of progenitors. 

5.4.2. Growth of the main HII region 

We aim with this section at understanding how the domi- 
nant region leads to a loss of information about the expan- 
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(d) Boosted Star 50 (e) Star 50 (f) Halo 50 



Fig. 14. Distribution of the HII regions sizes for the regions that merge with the dominant region as a function of redshift 
for the three kinds of ionizing sources and for both box sizes of 200 and 50 Mpc/h. Panels (a), (b) and (c) respectively 
represents the distribution for the Boosted Star, the Star and the Halo model for the 200 Mpc/h box, while the panels 
(d), (e) and (f) are for the same models but for the 50 Mpc/h box. The color code in the distributions is made so that the 
brightest red cell contains only a single HII region. Other red tones up to the brightest green designate cells populated by 
a couple HII regions up to a value corresponding to half the maximal value of the distribution. Finally the purple tones 
(from the darkest to the brightest) encode the largest values of the distribution. The evolution of the average radius 
for the regions that merge with the dominant region and the evolution of the radius of the main region are represented 
respectively with the solid and dashed white lines. 



sion process of others regions. By looking at the evolution 
of the properties of the regions that merge with the dom- 
inant one, we will evaluate how each model is resistant to 
the emergence of this main HII region. 



Fig. 14 presents the evolution of the radius distribution 
for HII regions that merge with the radius- dominant one 
for all models and both 200 and 50 Mpc/h boxes. With 
the white dashed curve we show again the evolution of the 
radius of the dominant region and with the solid white line 
we show the evolution of the averaged radius of the HII 
regions that merge with it only. The color coding of the 
distributions is identical to the one used in the figures pT 
and[l3l 

Immediately we see that, from a certain redshift, every 
model presents a radius distribution which is representative 
of the whole distribution of all HII regions, as discussed 
previously and shown in figure 11 In other words, from a 
certain moment onwards, the major part of the present HII 
regions are regions that will merge with the main region in 
the next snapshot. Thus this moment can be seen as the 
time when the main region imposes its domination. Then 
the later this moment appears, the longer the individual 
expansion process of ionized region can be tracked. In all 



models and at both resolution this moment is broadly co- 
incident with zbkg- from this redshift any region (in terms 
of size) can be incorporated into the background. 

It is also noteworthy to see how this main region is built 
before it begins to impose its domination. For instance both 
H50 and H200 models, the earliest stage of the buildup does 
not include any merger and the growth is purely driven by 
the inner source until z=9.5 (resp. 11.8) for the H200 (resp. 
H50) model. The equivalent region in the S200 model is 
detected later, with a large radius (~ 3 Mpc/h instead of 
0.5 Mpc/h for H200) and incorporates immediately other 
regions with similar sizes. At higher resolution, the radius 
dominant region of S50 is always larger than its H50 coun- 
terpart and also incorporate other relatively similar region 
quite early (at z ^ 14 for instance). Clearly the stronger 
emissivity of individual stellar source produce larger pro- 
genitors for the S models than for the halo-based ones, pro- 
moting a merger-driven growth sooner. Meanwhile the SB 
models exhibits the expected and quite different evolutions 
for the dominant region: the large emissivity correction at 
early times creates a large dominant region at the earliest 
stages and it quickly incorporates other regions also 'puffed 
up' by the boost and therefore belonging to the same class 
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of large radius. In particular, SB50 present a whole succes- 
sion of large regions, close to the white evolution track of 
the dominant one. For the latter, it implies that its radius is 
always larger than found in the S50 and H50 models and is 
driven by mergers much sooner. Previously, we found that 
high resolution tend to diminish the discrepancies of the SB 
model but in the current case, the higher level of cluster- 
ing in the 50 Mpc/h box leads to a buildup scenario of the 
main region that is arguably more even more distinguish- 
able compared to H and S models than in the 200 Mpc/h 
experiments. 

6. Summary and conclusion 

We developed a new methodology based on merger trees 
of HII regions to study the history of reionizations in 
cosmological simulations with full radiative transfer post- 
processed by the GPU-driven code AXON on RAMSES 
hydro dynamical snapshots. We applied systematically the 
technique in two sets of 200 Mpc/h and 50 Mpc/h simula- 
tions, where each set involved three kinds of source model: 

• a halo model where halos act as sources of photons with 
an emissivity proportional to their mass, 

• a stellar model where 'star' particles produced by the 
cosmological code act as the sources, with an emissiv- 
ity corrected to complete the reionization at the same 
redshift as the halo model, 

• a boosted stellar model, where the same 'star' parti- 
cles are used but with time-dependent and decreasing 
correction on the emissivity to reproduce the converged 
emission of UV photons at each instant. 

Corrections are necessary as the self-consistent production 
of 'star particles' by the hydrodynamical code is resolution- 
dependent and not fully converged at these resolutions. The 
emissivities were tuned to produce a similar ionization frac- 
tion evolution at both resolutions. HII regions were de- 
tected thanks to a FOF algorithm and linked along time 
within a merger tree structure. 

Firstly, the three models present an evolution compara- 
ble in terms of global features. Indeed, the evolution of the 
optical depth and ionized fraction present similar shapes in 
all experiments. 

The evolution of the number of present HII region 
reaches a maximum at the same similar redshift for all mod- 
els and for both box sizes (z=8.5 and z=9 for the 200 and 
50 Mpc/h experiments). However the study of the separate 
evolution of the different kinds of HII regions shows that 
the intense episodes of mergers occur earlier in the boosted 
Star model while it occurs later and at the same moment 
in the Star and Halo models in the 200 Mpc/h box. Also 
the boosted Star model is sensitive to recombination, due 
to sources unable to sustain the existing HII regions while 
only few recombinations occur in the Star model and no 
one in the Halo model in the 200 Mpc/h. In the 50 Mpc/h 
simulations, episodes of mergers are similarly distributed 
in the three models, whereas increased densities are likely 
to produce the larger population of recombining patches in 
the Star experiments. 

We also investigated the evolution with redshift of the 
size distribution of HII regions. All models present the 
emergence of a main region in size from a certain redshift 
onwards. We also have found that the evolution of the radii 



distribution is typically related to the evolution of the ion- 
izing source emissivities in each model: a large dispersion of 
possible radii with small structures in the halo-based mod- 
els, a similar distribution for the stellar model with a lack 
of small regions due to stronger and sparser emitters and a 
distribution strongly skewed toward large regions at early 
times in the boosted stellar experiments. The discrepancies 
decrease at higher resolution but remain, especially for the 
last model. We have then be able to show that this radius- 
redshift correlation, imprinted in the distribution of new 
regions is somehow memorized and kept in the distribution 
of growing and merging ionized patches. 

The evolution of the number of parents involved in the 
percolation of HII regions has also been studied. In all 
experiments, merger-dominant (with 10 or more progen- 
itors) regions exist and among them one corresponds to 
the radius-dominant one and drives the reionization. Its 
domination is delayed in both halo and stellar model com- 
pared to the boosted one. As an illustration, we investigated 
the size distribution of HII regions assimilated by the main 
one and found that indeed it becomes representative of the 
whole distribution first in the boosted stellar model, than 
in the star model and finally in the halo one. Hence individ- 
ual reionizations can be tracked on a longer period in the 
two latter models. Increasing the resolution reduces again 
the differences, thanks to an increased convergence between 
the three types of sources. 

Globally, we can say that the Star and the Halo models 
have similar histories even though the Star model lack small 
scale power and therefore local merger events. Meanwhile 
the boosted Star model presents its own reionization his- 
tory, where we see the early emergence of a dominant HII 
region in size that concentrates rapidly the merger pro- 
cess. On the other hand, both other models show individual 
growths for HII regions that can be tracked individually on 
a longer time and resist to the domination of the main re- 
gion later on during the reionization. In the 50 Mpc/h box, 
when the spatial resolution is increased, the histories of the 
three models become comparable even if the boosted Star 
model still presents the domination of an early main region 
that concentrates the merger process. 

Clearly, the large scale experiments show that the lack 
of convergence in the self- consistent formation of star parti- 
cles leads to substantial differences in the HII regions prop- 
erties and their evolution during the percolation process, 
even though they present satisfying properties regarding 
their global 'integrated' features. At higher resolution (cor- 
responding to 50 Mpc/h boxes in our cases), we find that 
overall the discrepancies are reduced and can be seen as the 
result of smaller corrections applied to sources and a greater 
match between the halo population and the 'stellar' one. A 
question also remains on the quantitative impact of the rare 
'emitters' more likely in the large boxes or the lack of den- 
sity field variance in the 50 Mpc/h. Future experiments will 
allow to assert them even though we can predict e.g. that 
rare emitters should be relevant for the rise of the domi- 
nant regions whereas density field variance on small scales 
should impact the rate of local mergers. 

We aim at systematically applying this technique in the 
future. For instance, it provides a reproducible analyzes 
framework to compare simulations, numerical techniques 
and the impact of physical ingredients on the morpho- 
chronology of the reionization. Furthermore, it puts the 
emphasis on individual HII regions, leading to an analyzes 
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in terms of reionization-s instead of The Reionization. In 
the context of galaxy formation, it may provide insights on 
the proximity effects which impacts the properties of the 
reionization as bei ng seen by one or a group of galaxy (see 
Ocvirk et ar]|2Q12| submitted) . 
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Appendix A: FOF algorithm 

A.l. Implementation 

Our FOF algorithm proceeds by scanning all cells from the 
box and testing at each iteration the ionization status of 
the cell being explored. The algorithm simply proceeds as 
follow: 
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Fig. A.l. Illustration of the friend- of- friend algorithm (see 
Section [a] for details). White and black cells correspond to 
ionized and neutral cells respectively. 



• If an ionized cell is encountered, it is given an ID cor- 
responding to the ionized region being tested and the 
exploration of its neighbors begins : 

1. If a neighbor is not ionized, the cell is marked as 
visited. The parent cell is then reconsidered and the 
ionization status of the other near neighbors not al- 
ready marked is tested. 

2. If the neighbor is ionized, it is given an ID corre- 
sponding to the ionized region being tested and we 
mark the cell as visited. The exploration of its own 
neighbors begins. The HII regions IDs diffuse from 
near neighbor to near neighbor. 

• When all the neighbors of an ionized cell have been ex- 
plored, we go back to its parent cell and continue the 
exploration of its near neighbor which were not previ- 
ously visited. The identification stops when the starting 
cell is visited again and all its near neighbors are marked 
as visited. 

The exploration of the cosmological box proceeds un- 
til another starting ionized cell is found and which is not 
already marked as visited. At this moment, a new identifi- 
cation starts with a new identification number. 



In Figure |A.1[ an illustration of the algorithm in 2 di- 
mensions is presented. In this figure, black cells correspond 
to neutral cell and white cells to ionized cells with a ion- 
ized fraction x > 0.5. In this example, 6 cells belongs to the 
same HII region, with numbers 17, 18, 25, 32 and 33 and 
all of them should receive the same identification number. 
Following the algorithm described above, the cells would be 
visited in the following order : 17-18-19-18-25-32-33-32-25- 
18-17. At this stage, the algorithm recover the starting cell 
with all its ionized neighbors marked as visited, it stops. 

To implement this method, a recursive scheme could be 
used as well as linked lists. We chose the second option 
as it limits the use of large stacks of memory. Let us also 
mention that for boxes of 1024^, we have parallelized the 
identification on the large box by cutting it in 16 sub-boxes. 
We thus do the identification on this 16 sub-boxes in par- 
allel to accelerate the calculation. We reassemble the entire 
box after and merge the IDs. 
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Fig. A. 2. Representation of the radius distribution of all HII regions at three different redshifts for the three simulations 
for the 50 Mpc/h box. We compare the difference in the radius distribution according to the ionization threshold for the 
detection of the HII region with the FOF procedure. We represent the volume of one single cell with the vertical red line. 



A.2. Effect of the neutral ionization threshold 

The FOF algorithm assumes an ionization threshold in or- 
der to decide if a cell is ionized or not. As often used in the 
literature, we have chosen a value of x > 0.5 to consider a 
cell as ionized. In order to test the impact of the threshold's 
value for the identification of the HII regions we have also 
tried other values for the ionization fraction in the FOF 
procedure. 



The figure A. 2 represents the radius distribution of the 
HII regions for three different instants. The distributions 
are done for the three model of ionizing sources for the 50 
Mpc/h box size. We compare the radius distributions ac- 
cording to three different ionizing thresholds of x > 0.3, 
X > 0.5 and x > 0.9 to assume that a cell is ionized. The 



three distributions are very similar for all models and at 
each instant regardless of the chosen ionizing threshold. The 
most significant variation stands for small regions with ra- 
dius below 2 X 10~^ Mpc/h that are detected in greater 
number with x > 0.9. Meanwhile, the distributions with 
X > 0.3 and x > 0.5 are almost superimposed in every 
cases. Additionally, we used the mer ger t ree properties to 
plot the same distributions in figure |A.3| but dissociating 
the contributions of each type of regions. Thus new regions 
are in blue, expanding regions in red, regions resulting from 
mergers in green and the regions that will recombine in yel- 
low. This time, we have just represented the distribution for 
both X > 0.5 and x > 0.9 ionization threshold, that previ- 
ously exhibited the greater difference. We see that globally 



19 



Jonathan Chardin et al.: Merger trees of HII regions 
Boosted Star Star 



Halo 






Fig. A. 3. Radius distribution of all HII regions at three different redshifts for the three simulations for the 50 Mpc/h box. 
This time, we dissociate the distribution in the different counterpart of the type of HII region considered. We compare 
the difference in the radius distribution according to both ionization threshold of x > 0.5 and x > 0.9 for the detection 
of the HII region with the FOF procedure. The color codes are as follows: Blue: the new regions, Red: the expanding 
regions. Green: the regions resulting from merger and Yellow: the regions that will recombine. In addition we represent 
the volume of one single cell with the vertical black line. 



that the families remain unchanged and therefore that the 
merger tree properties are conserved even when varying the 
ionization threshold. The only difference being that with 
X > 0.9 the FOF procedure tends to detect more small re- 
combining regions than with x > 0.5 in the boosted star 
and star models. 

Therefore the only remarkable effects by varying the 
threshold is for high threshold values within recombining 
regions, that represent a tiny number and volume fraction 
of the HII regions network. It does not come as a surprise 
though, since recombination occur in already ionized re- 



gions and any fluctuations within, combined with a high de- 
tection threshold lead to a subdivision of regions in smaller 
ones. One might even consider that such effect should be 
avoided since it increases branches in the merging tree with 
regions that may be short lived or unstable. As a conclu- 
sion, we can thus consider that the value of x > 0.5 used in 
the paper for the identification of the HII regions is a good 
compromise, in addition of being the one commonly used 
in the literature. 
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t t+i 

Fig. B.l. Merger tree algorithm. All neutral cells are black with an identification number of -1. All ionized cells are 
white with an identification number. At time t, we have the ionized region identified with an identification number of 1. 
We look forward at time t+1 where cells of this ionized region are located (shown in red). Here ID 1 at time t is linked 
to ID 2 at time t+1. 



Appendix B: Merger tree algorithm 

The implementation of the mer ger t ree algorithm is per- 
formed as follow (see also Figure B.l): 



1. 



2. 



3. 



We look at an instant t of the simulation where are 
located the cells of a particular ionized region. They all 
share the same ID. 

Then, we extract the same cells at time t-\-l. These cells 
may have different IDs. 

The most common ID at time t + 1 is kept and is linked 
to the ID investigated at t. 



This scheme is repeated for all HII regions between the 
snapshot t and t We finally reproduce this process be- 
tween each snapshots. 
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